---
title: "Experiment 4"
output: word_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r load libraries}
library(readr)
library(reshape2)
library(plyr)
library(ggplot2)
library(mediation)
library(tidyverse)
```

```{r import}
fulldata <- read_csv("~/Desktop/Studies/Force/analyses/FINAL/Exp.4/exp4_foranalysis.csv", 
    col_names = FALSE)


```
##Clean data
```{r clean_data echo=False}
#Name the columns
names(fulldata)<-c("sub","duration","scenario","cong_incong","force","norm","value")

#Remove rows with any missing data
cleandata<-fulldata[complete.cases(fulldata), ]

#Recode scenario to be 0, 1, 2
cleandata <- cleandata %>%
    mutate(scenario = scenario - 1)

#Specify that the variables are factors
cleandata$scenario<-factor(cleandata$scenario,labels=c("lab","taxi","track"))
cleandata$sub<-factor(cleandata$sub)
cleandata$cong_incong<-factor(cleandata$cong_incong, labels=c("cong","incong"))

#separate out scenarios for other analyses
labdata=cleandata[which(cleandata$scenario=="lab"),]
taxidata=cleandata[which(cleandata$scenario=="taxi"),]
trackdata=cleandata[which(cleandata$scenario=="track"),]

```
##Descriptives
```{r descriptives}
mosaic::favstats(~force, groups=cong_incong, data=cleandata)

#get descriptives
force_desc <- ddply(cleandata, c("scenario","cong_incong"), summarise,
               N    = length(force),
               mean = mean(force),
               sd   = sd(force),
               se   = sd / sqrt(N)
)
print(force_desc)

norm_desc <- ddply(cleandata, c("scenario","cong_incong"), summarise,
               N    = length(norm),
               mean = mean(norm),
               sd   = sd(norm),
               se   = sd / sqrt(N)
)
print(norm_desc)

value_desc <- ddply(cleandata, c("scenario","cong_incong"), summarise,
               N    = length(value),
               mean = mean(value),
               sd   = sd(value),
               se   = sd / sqrt(N)
)
print(value_desc)

```
##Final Plots

```{r final visualizations for publication}

###FINAL BAR CHART####
#categorize study 2 data
exp4_val_summary<-mosaic::favstats(~value,data=cleandata)
exp4_norm_summary<-mosaic::favstats(~norm,data=cleandata)

cleandata$val_cat<-ifelse(cleandata$value >= exp4_val_summary$Q3,cleandata$val_cat<-"High Value",ifelse(cleandata$value <= exp4_val_summary$Q1, cleandata$val_cat<-"Low Value",cleandata$val_cat<-"Medium Value"))

cleandata$norm_cat<-ifelse(cleandata$norm >= exp4_norm_summary$Q3,cleandata$norm_cat<-"High Unusualness",ifelse(cleandata$norm <= exp4_norm_summary$Q1, cleandata$norm_cat<-"Low Unusualness",cleandata$norm_cat<-"Medium Unusualness"))


#Create factor out of bin factors and reorder
cleandata$val_cat<-factor(cleandata$val_cat, ordered=TRUE, levels=c("Low Value", "Medium Value", "High Value"))
cleandata$norm_cat<-factor(cleandata$norm_cat, ordered=TRUE, levels=c("Low Unusualness", "Medium Unusualness", "High Unusualness"))

#Create condition factor for graph labels
cleandata$Condition<-dplyr::recode(cleandata$cong_incong, incong="Normal Alternative", cong="Abnormal Alternative")
cleandata$Condition<-factor(cleandata$Condition, levels=c("Normal Alternative", "Abnormal Alternative"))

####Create final bar chart####
cond_legend_title<-"Alternative Action"
exp4_fig2a<-cleandata %>%
  ggplot(aes(Condition, force, fill=Condition)) +
  geom_bar(stat="summary", fun.y='mean', position='dodge') +
  geom_errorbar(stat='summary', width=.2,
                 position=position_dodge(.9)) +
  geom_point(aes(Condition, color=Condition), size=2,
             position=position_jitterdodge(dodge.width=.9,jitter.width = .8), alpha = .5) +
             scale_color_manual(cond_legend_title, values = c("Normal Alternative" = "#8a433e", "Abnormal Alternative" = "#017073")) +
  scale_fill_manual(cond_legend_title, values = c("Normal Alternative" = "#F8766D", "Abnormal Alternative" = "#00BFC4")) +
  xlab("") +
  ylab(str_wrap("How strongly do you agree that the agent was forced?", 75)) +
  scale_y_continuous(breaks=waiver(), labels = str_wrap(c("Strongly Disagree","25","50","75","Strongly Agree"),10)) + 
  theme_bw() +
  theme(strip.text = element_text(size=14),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title = element_text(size=20, vjust=0.35),
        panel.grid = element_blank(),
        legend.position = "none")

ggsave("exp4_fig2a.jpg",height = 10, width= 7, dpi=600)

###BINED BAR CHARTS####
val_legend_title<-"Value Ratings"
unus_legend_title<-"Unusualness Ratings"
exp4_fig2b<-cleandata %>%
  ggplot(aes(val_cat,force, fill=val_cat)) +
  geom_bar(stat="summary", fun.y='mean', position='dodge') +
  geom_errorbar(stat='summary', width=.2,
                 position=position_dodge(.9)) +
  geom_point(aes(val_cat, color=val_cat), size=2,
             position=position_jitterdodge(dodge.width=.9,jitter.width = .8), alpha = .5) +
  scale_color_manual(val_legend_title,values = c("Low Value" = "#919191", "Medium Value" = "#6e6e6e", "High Value" = "#505050")) +
  scale_fill_manual(val_legend_title,values=c("Low Value" = "#c3c3c3", "Medium Value" = "#A9A9A9", "High Value" = "#808080")) +
  xlab("") +
  ylab(str_wrap("How strongly do you agree that the agent was forced?", 30)) +
  scale_y_continuous(breaks=waiver(), labels = str_wrap(c("Strongly Disagree","25","50","75","Strongly Agree"),10)) + 
  theme_bw() +
  theme(strip.text = element_text(size=14),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title = element_text(size=20, vjust=0.35),
        panel.grid = element_blank(),
        legend.position = "none")

ggsave("exp4_fig2b.jpg",height = 5, width= 7, dpi=600)

exp4_fig2c<-cleandata %>%
  ggplot(aes(norm_cat,force, fill=norm_cat)) +
  geom_bar(stat="summary", fun.y='mean', position='dodge') +
  geom_errorbar(stat='summary', width=.2,
                 position=position_dodge(.9)) +
  geom_point(aes(norm_cat, color=norm_cat), size=2,
             position=position_jitterdodge(dodge.width=.9,jitter.width = .8), alpha = .5) +
             scale_color_manual(unus_legend_title,values = c("Low Unusualness" = "#919191", "Medium Unusualness" = "#6e6e6e", "High Unusualness" = "#505050")) +
  scale_fill_manual(unus_legend_title,values=c("Low Unusualness" = "#c3c3c3", "Medium Unusualness" = "#A9A9A9", "High Unusualness" = "#808080")) +
  xlab("") +
  ylab(str_wrap("How strongly do you agree that the agent was forced?", 30)) +
  scale_y_continuous(breaks=waiver(), labels = str_wrap(c("Strongly Disagree","25","50","75","Strongly Agree"),10)) + 
  theme_bw() +
  theme(strip.text = element_text(size=14),
        axis.text.x = element_text(size=14),
        axis.text.y = element_text(size=14),
        axis.title = element_text(size=20, vjust=0.35),
        panel.grid = element_blank(),
        legend.position = "none")

ggsave("exp4_fig2c.jpg",height = 5, width= 7.7, dpi=600)
```

#Final analyses
```{r testing effect of condition echo=FALSE}
#Predict force judgments by condition, with a random intercept for scenario.
fullmod<-lme4::lmer(force ~ cong_incong + (1|scenario), data=cleandata)
summary(fullmod)

#New model excluding condition
modnocond<-lme4::lmer(force ~ (1|scenario), data=cleandata)
summary(modnocond)

#Compare the models
anova(fullmod,modnocond)

```

```{r subject level multiple regression echo=FALSE}
#Predict force judgments with unusualness (norm) and value judgments, with a random intercept for scenario. 
fullmod<-lme4::lmer(force ~ norm + value + (1|scenario), data=cleandata)
summary(fullmod)

#New model excluding unusualness judgments
modnonorm<-lme4::lmer(force ~ value + (1|scenario), data=cleandata)
summary(modnonorm)

#Compare to original model
anova(fullmod,modnonorm)

#New model excluding value judgments
modnoval<-lme4::lmer(force ~ norm + (1|scenario), data=cleandata)
summary(modnonorm)

#Compare to original model
anova(fullmod,modnoval)

#Rerun original model with condition included
fullmod_withcond<-lme4::lmer(force ~ norm + value + cong_incong + (1|scenario), data=cleandata)
summary(fullmod_withcond)

#Compare to original model
anova(fullmod,fullmod_withcond)

```


